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Abstract 

A low-dimensional multidimensional scaling example is used to illustrate properties 
of the stress loss function and of different iteration methods. 
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Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome, and some would be really benificial. For example, I only use 
base R graphics, nothing more fancy, because base is all I know. The directory deleeuw- 
pdx.net/pubfolders/twoPoints has a pdf version, the complete Rmd file with all code chunks, 
the bib hie, and the R source code. 
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1 Theory 

In Multidimensional Scaling (MDS) we minimize a loss function, called Kruskal’s stress, which 
is defined by (Kruskal (1964a), Kruskal (1964b)) 


a(Z) := 


inn 

2 

Z i =1 3 = 1 


d zj (Z)f 


( 1 ) 


over the n x p configurations Z. Here the weights W and the dissimilarities A are known 
symmetric, non-negative, and hollow matrices of numbers. The matrix D{Z) in (1) has the 
(Euclidean) distances between the n points in the configuration. Thus if e t are unit vectors, 
with zeroes everywhere, except for element i, which is one, and Aij := (e* — ej)(e; — ef )', then 
djjiZ) = tr ZAijZ'. 

Stress is a complicated function with potentially a large number of local minima and saddle 
points. In order to study the behavior of stress we shall look at configurations of the form 
Z = aX + /3Y, with X and Y known, fixed, and centered configurations. This makes stress 
a function of the two variables a and (d, and we can use standard contour and perspective 
plots to study the function. This extends earlier work by De Leeuw (1993). We assume linear 
independence, in the sense that aX + /3Y = 0 if and only if a — /3 — 0. 

First some simplications. Define the 2x2 matrices 


Va := 


tr XAijX' tr XA t] Y' 
tr YAijX' tr YA l3 Y' ’ 


and define 7 as the vector with elements a and fd. Now 


a(y) 


n n _ 

Z E w iAj\fr%Yr+ ~VK*7, 


i =1 3 = 1 


where we have assumed for convenience that 


1 

2 


E E w ii s ij = *• 


i=i j =1 


and where 

n n 

i= 1 j =1 

Note that if all w %3 are one, then 


K* 


2 n 


X'X 

Y'X 


X'Y 

Y'Y 


( 2 ) 


We now make a change of variables, using the Cholesky decomposition W* = S'S, with S 
upper-triangular. Define 6 := S 7 and U l3 := [S')^ l VijS~ l . Then 
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( 3 ) 


n n _ i 

v{°) = l ~Y,Yl WijSijy/O'Uijd + -6'6. 

i= i j =i 

In MDS we usually (De Leeuw (1977)) also define the matrix-valued function 


5. 


i= 1 j =1 d l j(9) 


(4) 


with dij{9 ) := sJd’UijO and the function p($) := 9'B(9)9. Definition (4) can be extended 
using subgradients if dij(9) = 0 for some i, j (see De Leeuw (1977)). For our purposes we can 
define B{9) more generally by omitting the terms with dij(9) = 0. The definition of p allows 
us to write 


a(9) = 1 - p(6) + X -9'9. (5) 

If dijiQ) > 0 for all i , j such that WijSij > 0 then stress is differentiable at 9 and T>a(9) = 
9 — B(9)9, so that stationary points are defined as solutions of 9 = B{9)9. The more general 
definition is 9 G dp{9). De Leeuw (1984) shows that stress is always differentiable at local 
minima. 

At a stationary point 9 is an eigenvector of B{9 ) with eigenvalue equal to one. If this is the 
largest eigenvalue, then 9 actually gives the global minimum of stress (De Leeuw, Groenen, 
and Mair (2016)). Also observe that at a stationary point p(9) = 9’9, and thus (5) implies 
9'9 < 2. All stationary values are within a circle or sphere with radius \/2. 

The smacof algorithm (De Leeuw (1977), De Leeuw and Mair (2009)) is the iterative algorithm 


e (k+l) = £(0(fc))0(fc)_ 


( 6 ) 


Note that the algorithm is self-scaling, in the sense that all points on a ray through the origin 
have the same update. Thus B(X9)(X9) = B(9)9 for all A ^ 0. 

If stress is differentiable at 9 then V 2 a{9) = I — H(9), where 


me) 


n n X 

Oij 


W. 


*3 


i =1 3 =1 


dij{9) 


U^Un \ 
9'Uij9 / 


(7) 


Note that Vp{9) = B{9)9 and V 2 p{9) = H{9). Also note that H(X9) = |A| l H{9 ), and thus 
for any 9 there is a A (9) such that T> 2 cr(X9) is positive definite for all A > A (9). 

The smacof algorithm convergences to a solution 9 of the stationary equations, with a linear 
convergence rate equal to the largest eigenvalue of H(9) (see De Leeuw (1988)). 

Note that V 2 a{9)9 = 9 , which means that Newton’s method takes the form 
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9 (k+i) = (/ _ H(6 {k) y l B(6 {k) )6 {k \ 


( 8 ) 


Note that Newton is definitely not self-scaling. Equation ( 8 ) also suggest a natural safeguarded 
version of Newton’s method, using / — pH, where 0 < 77 < 1. For 77 = 0 this is a smacof step, 
for 77 — 1 this is a Newton step. 

On the line through 9 and the origin the function a is a convex quadratic. Thus at every 
point, except at the origin, there is at least one direction of descent, and consequently stress 
has only a single local maximum at 9 = 0, equal to one. In the differentiable case this is also 
clear from T> 2 a{9)9 = 9, which says that the Hessian has at least one positive eigenvalue. 

Because stress is an even function, with a(9) = —a(9) for all 9 , minima come in pairs. 


2 Example 

Our first example has n = 4, all weights equal to one, and all dissimilarities equal. The 
same example has been analyzed by De Leeuw (1988), De Leeuw (1993), Trosset and Mathar 
(1997). For this example the global minimum in two dimensions has its four points in the 
corners of a square. That is our X, which has stress 0.02859548. Our Y is another stationary 
point, which has three points in the corners of an equilateral triangle and the fourth point 
in the center of the triangle. Its stress is 0.0669873. Another way of looking at the two 
configurations is that X are four points equally spaced on a circle, and Y are three points 
equally spaced on a circle with the fourth point in the center of the circle. De Leeuw (1988) 
erroneously claims that Y is a non-isolated local minimum of stress, but Trosset and Mathar 
(1997) have shown there exists a descent direction at Y, and thus Y is actually a saddle point. 
Of course the stationary points defined by A" and Y are far from unique, because we can 
distribute the four points over the various corners in many ways. 

The example is chosen in such a way that there are non-zero a and [3 such that dui^aX + /3Y) = 
0. In fact di 2 is the only distance that can be made zero by a non-trivial linear combination. 

Note that we have used o for three different functions. The first one with argument Z is 
defined on configuration space, the second one with argument 7 on coefficient space, and the 
third one with argument 9 also on coefficient space. This is a slight abuse of notation, rather 
innocuous, but we have to keep it in mind. 

From lemma lwe see that Va(X) = Va(Y) = 0 then Da(l,0) = T>a(0, 1) = 0. Thus 
stationary points in configuration space are preserved as stationary points in coefficient space, 
but the reverse implication may not be true. If T> 2 <j(X) and V 2 a(Y) are positive semi-definite, 
then so are T> 2 (j(l,0) and T> 2 cr( 0,1). Thus local minima are preserved. But it is entirely 
possible that T> 2 a(X) and/or T> 2 afiY) are indefinite, and that T> 2 a( 1,0) and/or T> 2 a( 0,1) are 
positive semi-definite. Thus saddle points in configuration space can be mapped into local 
minima in coefficient space. As we will see this actually happens with Y, the equilateral 
triangle with center, in our example. 
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2.1 Global Pictures 

We first make a global perspective plot, over the range (—2.5, +2.5). 



theta 1 


Figure 1: Global Perspective 

We see the symmetry, following from the fact that stress is even. We also see the local 
maximum at the origin, where stress is not differentiable. Also note the ridge, where d 12 ($) = 0 
and where stress is not differentiable either. The ridge shows nicely that on rays emanating 
from the origin stress is a convex quadratic. Also, far away from the origin, stress globally 
behaves very much like a convex quadratic (except for the ridge). Clearly local minima must 
be found in the valleys surrounding the small mountain at the origin, all within the sphere 
with radius y/2. 

Figure 2 is a countour plot of stress over (—2, +2) ® (—2, +2). The red line is {0 | di2(0) = 0}. 
The blue line has the minimum of the convex quadratic on each of the rays through the origin. 
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Thus all local minima, and in fact all stationary points, are on the bine line. Note that the 
plot uses 9 to define the coordinate axes, not 7 = (a, (3). Thus there are no stationary points 
at (0, f) and (1,0), but at the corresponding points (1.3938468501, 0) and (1.040640449, 
0.8849253413) in the 6 coordinates (and, of course, at their mirror images). 

Besides the single local maximum at the origin, it turns out that in this example there are at 
least five pairs of stationary points. Or, more precisely, 1 have not been able to find more 
than five. Each stationary point 6 has a mirror image —6. Three of the five are local minima, 
two are saddle points. Local minima are plotted as blue points, saddle points as red points. 



Figure 2: Global Contour 
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2.2 First Minimum 


We zoom in on the first local minimum at (1.0406404488, 0.8849253415). Its stress is 
0.0669872981 and the corresponding configuration is in figure 3. 
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dim 1 


Figure 3: Configuration First Minimum 

Note that this local minimum corresponds with the equilateral triangle with center, which is 
a saddle point in configuration space (Trosset and Mathar (1997)). The eigenvalues of B{9) 
are (1.3686346146, 1) and those of the Hessian / — H{6) are (1, 0.0817218104). The area of 
the contour plot around the stationary value is in figure 4. 
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Figure 4: Contour Plot First Minimum 
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2.3 Second Minimum 

The second local minimum (which is the global minimum) at (1.3938468501, 0) has stress 
0.0285954792 and the corresponding configuration is in figure 5. 
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Figure 5: Configuration Second Minimum 

The eigenvalues of B(6) are (1.1362798549, 1) and those of the Hessian / — H{9 ) are (1, 
0.3743105348). The area of the contour plot around the stationary value is in figure 6. 
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Figure 6: Contour Plot Second Minimum 

2.4 Third Minimum 

The third local minimum at (0.1096253665, 1.3291941942) has stress 0.1106125366 and the 
corresponding configuration is in figure 7. 
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Figure 7: Configuration Third Minimum 

The eigenvalues of B(9) are (1.5279386667, 1) and those of the Hessian / — H{9 ) are (1, 
0.2362078313). The area of the contour plot around the stationary value is in figure 8. 
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Figure 8: Contour Plot Third Minimum 

2.5 First Saddle Point 

The saddle point at (0.3253284471, 1.2916758288) has stress 0.1128674774 and the corre¬ 
sponding configuration is in figure 9. 
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Figure 9: Configuration First Saddlepoint 

The eigenvalues of B{9) are (1.7778549251, 1) and those of the Hessian / — H{6) are (1, 
-0.3110880057). The area of the contour plot around the stationary value is in figure 10. 
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Figure 10: Contour Plot First Saddlepoint 


2.6 Second Saddle Point 

The saddle point at (1.1238371005, 0.7762045565) has stress 0.067248329 and the correspond¬ 
ing configuration is in figure 11. 
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Figure 11: Configuration Second Saddlepoint 

The eigenvalues of B{6) are (1.4111961819, 1) and those of the Hessian / — H{6) are (1, 
-0.0841169222). The area of the contour plot around the stationary value is in figure 12. 
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Figure 12: Contour Plot Second Saddlepoint 

3 Regions of Attraction 

3.1 Smacof 

We use the smacof () function from the code in the appendix with 100 different starting 
points of 9 , equally spaced on the circle. Figure 13 is a histogram of the number of smacof 
iterations to convergence within le — 15. In all cases smcof converges to a local minimum in 
coefficient space, never to a saddle point. Figure 14 shows which local minima are reached 
from the different starting points. This shows, more or less contrary to what Trosset and 
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Frequency 


Mathar (1997) suggest, that non-global minima can indeed be points of attraction for smacof 
iterations. 

Histogram of iterations 
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Figure 13: Histogram Number of Smacof Iterations 
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Figure 14: Path Endpoints of Smacof Iterations 

3.2 Newton 

We repeat the same exercise with Newton’s method, which converges from all 100 starting 
points. In higher dimensions we may not be so lucky. The histogram of iteration counts is in 
figure 15. It shows in this example that smacof needs about 10 times the number of iterations 
that Newton needs. Because smacof iterations are much less expensve than Newton ones, 
this does not really say much about computing times. If we look at figure 16 we see the 
problem with non-safeguarded Newton. Although we have fast convergence from all 100 
starting points, Newton converges to a saddle point in 45 cases. 
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Figure 15: Histogram Number of Newton Iterations 
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Figure 16: Path Endpoints of Newton Iterations 

4 Another Look 

Remember that p{9) = 9'B{9)9. Thus a(X9) = 1 — A p{9) + 2 9'9 1 and 

Thus we can minimize a over 9 by maximizing p over the unit circle S := {9 \ 9'9 = 1}. 
This is a nice formulation, because p is norm, i.e. a homogeneous convex function of 9. 
Consequently we have transformed the problem from unconstrained minimization of the DC 
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function (i.e. difference of convex functions) stress to that of maximization of a ratio of norms. 
In turn this is equivalent to maximization of the convex function p over the unit circle, or, 
again equivalently, over the unit ball, a compact convex set. This transform was first used in 
MDS by De Leeuw (1977), partly because it made the theory developed by Robert (1967) 
available. 

The levels sets {9 \ p[9) = k] are the p-circles defined by the norm p. The corresponding 
p-balls {6 j p{9) < k} are closed and nested convex sets containing the origin. Thus we 
want to hnd the largest p-circle that still intersects S. The similarity with the geometry of 
eigenvalue problems is obvious. 

In our example we know that the global optimum of stress is at (1.3938468501, 0), and 
if we project that point on the circle it becomes (1, 0). The corresponding optimal p is 
1.3938468501. Figure 17 gives the contourplot for p, with the outer p-circle corresponding 
with the optimal value. The fact that the optimal value contour is disjoint from the interior of 
S is necessary and sufficient for global optimality (Diir, Horst, and Locatelli (1998)). Notice 
the sharp corners in the contour plot, showing the non-diffentiability of p at the points where 
h 12 (0) = 0. We could also look for the minimum of p on the unit circle, which means finding 
the largest p-circle that touches S on the inside. Inspecting figure 17 shows that this will be 
a point where p is not differentiable, i.e. a point with d^iO) = 0. This minimum p problem 
does not make much sense in the context of multidimensional scaling, however, and it not- 
related directly to the minimization of stress. 
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Figure 17: Contour Plot for Rho 


5 A Final Look 

Now that we know that the MDS problem is equivalent to maximizing p on the unit circle, 
we can use nonlinear coordinates { 61 , 62 ) = (sin £, cos £) to reduce the problem to a one¬ 
dimensional unconstrained one in, say, “circle space”. Thus, with the same abuse of notation 
as for stress, p(£) := p(sin£, cos£), and we have to maximize p over 0 < £ < tt. 

In figure 18 we have plotted p as a function of 77 . There are blue vertical lines at the three 
local minima in coefficient space, red vertical lines at the stationary points, and a green 
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vertical line where d 12 (£) = 0. Note that in circle space stress has both multiple local minima 
and multiple local maxima. 
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Figure 18: One-dimensional Rho 

From lemma 2 we see that the second derivative T> 2 p (£) is equal to tr H(£) — p(£). For the 
three local minima in coordinate space we find second derivatives -0.111634069, -0.5217315599, 
-0.3150320881 in circle space, i.e. they are properly converted to local maxima. The two 
stationary points in coordinate space have second derivatives 0.4143740166, 0.1148897778, 
and are turned into local minima. 

For more general cases, with a basis of n configurations, we know from Lyusternik and 
Schnirelmann (1934) that a continuously differentiable even function on the unit sphere in 
R n has at least n distinct pairs of stationary points. 
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6 Discussion 


Although our results are based on a single small example, there are some general conclusions 
we can draw. 

Stess has no local maxima, except one at the origin. It has saddle points and non-global 
minima. If we reparametrize configuration space by using linear combinations of a fixed 
number of configurations we can “upgrade” some of the stationary point to local minima. If 
we remove homogeneity from the problem by working on the unit ball that can upgrade even 
more stationary points. 

In comparing smacof with Newton we have found, not surprisingly, that Newton uses fewer 
iterations but often converges to saddle points. 

The finding that saddle points in configuration space can correspond with local minima in 
coefficient space has an interesting implication for the smacof algorithm. We know that the 
set of starting points from which smacof converges to a saddle point has measure zero, in 
other words unless you start in a saddle point, you will almost certainly converge to a local 
minimum (Lee et al. (2016)). 

To illustrate this we repeat a calculation done first in De Leeuw (1988). Suppose we start 
smacof iterations with the Y + .001 * E, where E are random standard normals, and Y is 
the triangle with center (which has loss function value 0.0669872981). smacof starts with a 
loss function value for the perturbed Y of 0.066987868. It does not drop below 0.066 until 
iteration 1346. But then it rather quickly converges in iteration 1388 to 0.0285954792, the 
loss function value of the global minimum, four points in the corners of a square. On the 
other hand, if we minimize <r( 7 ) over 7 we have a non-zero probability of converging to the 
local minimum in coefficient space at ( 0 , 1 ), i.e. to Y. 

This indicates, perhaps, that smacof must be used with a somewhat higher precision than 
the default, and that testing for second derivative information is always a good idea, no 
matter what space we are working in. 

Minimizing stress over a plane spanned by two configurations may seem somewhat artificial 
and limited. But think of the situation in which we use configurations X and Y = T>a(X). 
Or, equivalently, X and B(X)X. In that case minimizing over the plane means computing 
the optimal step size of a gradient step or the optimum over-relaxation of smacof iterations, 
a problem first addressed perhaps in De Leeuw and Heiser (1980). Or the case in which Y is 
the Newton step, and optimizing over the plane is a stabilized version of Newton’s method. 

7 Appendix: Two Lemmas 

Here we present two lemmas that describe the change of coordinates from configuration space 
to coefficient space and then to circle space. The proofs, which are just simple computations, 
are omitted. 

Lemma 1: [On the Line] Suppose X and Y are nxm matrices and / is a twice-differentiable 
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function on M nxm . Define g(a, (3) := f(aX + /3Y ). Then 


Vg(a,f3) 


tr X'Vf(aX + f3Y) 
tr Y'T>f(aX + (3Y) ’ 


and 

V 2 g(a,l3) = 


EL,E” 1 EL,EL, XiiXujsjgs^aX + fSY) EL,E™,EL,EL, 

(“A' + 0Y) EL, E ”, EL, EL, 


V n -c-m v^n v^m 9/ 2 

2^j=l Z^ 7 = l 2^fc=l 2^=1 yij^ki 


It is clear how this result generalizes to linear combinations of more than two matrices. 

Lemma 2: [On the Circle] If / is a twice-differentiable function of two variables and 

g(x) = f(sin(x),cos(x )) then 

Vg(x) = y'Vf(z ), 

with z := (sin(x), cos(x)) and y := (cos(x), —sin(x)), and 

V 2 g{x) = y'V 2 f(z)y - z'Vf(z). 


More generally, if a differentiable function / on the unit sphere has a stationary point at 9 
then 9’9 = 1 and there is a multiplier A such that T>f{9) = A 6 . If / is homogeneous of degree 
one then 9"Df{9 ) = f{9 ) and thus at a stationary point A = f(9). If the function is twice 
differentiable and the stationary point is a local maximum then in addition v'V 2 f(9)v < A 
for all v’v = 1 and v'Q = 0. For a norm, i.e. a homogeneous convex /, this says equivalently 
that the largest eigenvalue of V 2 f{9) is less than or equal to f(9). Note that the second order 
necessary conditions for a local maximum become sufficient if the inequality is strict. 


8 Appendix: Code 

bmat <- function (a, b, x, y, delta) { 
bm <- matrix (0, 2, 2) 
hm <- matrix (0, 2, 2) 
z <- c(a, b) 
for (i in 1:4) { 
for (j in 1:4) { 
if (i == j) next 
uij <- uu (i, j , x, y) 
uz <- drop (uij 7 0 *°/ 0 z) 
dij <- sqrt (sum (uij * outer (z, z))) 
bm <- bm + (delta[i,j] / dij) * uij 

hm <- hm + (delta[i,j] / dij) * (uij - outer (uz, uz) / sum (z * uz)) 

> 

> 

return (list (b = bm, h = hm)) 


(aX + f3Y 
(aX + (3Y 
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> 


stress <- function (a, b, x, y, delta) { 
z <- c (a, b) 

bm <- bmat (a, b, x, y, delta) $b 

return (1 + sum(z 2) / 2 - sum (z * bm °/ 0 *% z)) 


rho <- function (a, b, x, y, delta) { 
z <- c (a, b) 

bm <- bmat (a, b, x, y, delta) $b 
return (sum (z * bm %*% z)) 


vv <- function (i, j, x, y) { 
a <- matrix (0, 2, 2) 
a[l, 1] <- sum ((x [i, ]- x[j,]) ~ 2) 
a [2, 2] <- sum ((y [i, ]- y [j ,]) ~ 2) 

a[l, 2] <- a [2, 1] <- sum ((x[i, ]- x[j,]) * (y[i, ]- y[j, ])) 
return (a) 


uu <- function (i, j, x, y) { 
n <- nrow (x) 

asum <- 2 * n * matrix (c (sum(x ~ 2 ), sum (x * y), sum (x * y), sum (y ~ 2)), 2, 2) 
csum <- solve (chol (asum)) 

return (t(csum) °/ 0 *°/o vv (i, j, x, y) %*% csum) 

} 

smacof <- function (a, b, x, y, delta, eps = le-10, itmax = 1000, verbose = TRUE) { 
zold <- c(a,b) 

bold <- bmat (a, b, x, y, delta)$b 

fold <- 1 + sum(zold 2) / 2 - sum (zold * bold %*% zold) 
itel <- 1 
repeat { 

znew <- drop (bold zold) 

bhmt <- bmat (znew[l], znew [2] , x, y, delta) 

bnew <- bhmt$b 

fnew <- 1 + sum (znew 2) / 2 - sum (znew * bnew %*% znew) 
if (verbose) { 
cat ( 

formate (itel, width = 4, format = "d") , 
formate ( 
fold. 
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digits = 10, 
width = 13, 
format = "f" 

), 

formate ( 
fnew, 

digits = 10, 
width = 13, 
format = "f" 

), 

"\n" 

) 

> 

if ((itel == itmax) |I (fold - fnew) < eps) 
break () 

itel <- itel + 1 
fold <- fnew 
zold <- znew 
bold <- bnew 

} 

return (list (stress = fnew, theta = znew, itel= itel, b = bnew, g = znew - bnew 7o*7o 


newton <- function (a, b, x, y, delta, eps = le-10, itmax = 1000, verbose = TRUE) { 
zold <- c(a,b) 

bhmt <- bmat (a, b, x, y, delta) 

bold <- bhmt$b 

hold <- diag(2) - bhmt$h 

fold <- 1 + sum(zold 2) / 2 - sum (zold * bold %*% zold) 
itel <- 1 
repeat { 

znew <- drop (solve (hold, bold 7o*7» zold)) 
bhmt <- bmat (znew[l], znew[2], x, y, delta) 
bnew <- bhmt$b 

hnew <- diag(nrow(bnew) ) - bhmt$h 

fnew <- 1 + sum (znew 2) / 2 - sum (znew * bnew 7o*7o znew) 
if (verbose) { 
cat ( 

formate (itel, width = 4, format = "d") , 
formate ( 
fold, 

digits = 10, 
width = 13, 
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format 


), 

formate ( 
fnew, 

digits = 10, 
width = 13, 
format = "f" 

), 

"\n" 

) 

> 

if ((itel == itmax) |I abs (fold - fnew) < eps) 
break () 

itel <- itel + 1 
fold <- fnew 
zold <- znew 
bold <- bnew 
hold <- hnew 

} 

return (list (stress = fnew, theta = znew, itel = itel, b = bnew, g = znew - bnew %*% 

} 

mprint <- function (x, d = 2, w = 5) { 

print (noquote (formate (x, di = d, wi = w, fo = "f"))) 

} 
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